9. Elastic Deformation#
Elastic deformation is computed using the ContactMechanics library.
Under the assumption that the surfaces are linear-elastic and isotropic, the deformation is computed using a boundary element method. Here, we are just interested in the displacement of the surface and not of the deformation within the material. This renders the 3D problem into a 2D problem that can be efficiently solved using fundamental solutions (so-called Green’s functions).
Theory#
Green’s Function Approach#
Assuming that the elastic deformation of the surfaces responds on a shorter timescale than the fluid-dynamic system, we demand the stresses \(\sigma\) in the material to be in equilibrium (\(\frac{\partial \underline{\sigma}}{\partial t} = 0\)), leading to the Navier-Lame equation.
After putting in Hooke’s law relating stress to strain or displacement \(u\), respectively, we obtain a non-linear PDE in three dimensions. With some tricks, this system can be analytically solved for specific cases, see e.g. Johnson, Cambridge University Press, 2012. Contact Mechanics.
To obtain the fundamental solution (Green’s function) for our approach, we solve the problem of a normal point force \(F_z\) acting on a infinite surface in \(x\)- and \(y\)-direction (infinite half-space) without any tangential or body forces. With this solution we can then calculate the normal displacement:
This is the essence of this approach: Solution of the displacement \(u_z(x,y)\) is obtained by convolution of the Green’s function \(G(x,y)\) with the pressure field \(p(x,y)\) over the domain \(A_G\). Note that \(G(x,y)\) only needs to be computed once.
The Green’s function is different for periodic and non-periodic boundary conditions. See the following sections for details on
periodic solution
non-periodic solution
semi-periodic solution (periodicity in one direction, non-periodic in the other)
Periodic Solution#
For periodic boundary conditions, the system can be solved in Fourier space. The resulting Fourier space Green’s function takes the form:
with the effective Young’s modulus \(E_{eff}\) and the wavevector \(\vec{q}\). Note that we assume infinite thickness of the elastic half space. We set \(\tilde{G}(q_0) = 0\), resulting in the mean displacement being zero.
Non-Periodic Solution#
The non-periodic, real-space Green’s function takes the form
with a singular point at the origin. To manage this singularity, the problem can be regularized by discretizing the domain into rectangles and to assume constant pressure over these rectangles.
By that, we obtain the real space solution for displacement of a general point \((x, y)\) due to a uniform pressure acting on a rectangular area \(2a \cdot 2b\) with center in the origin
with the analytic solution
Semi-Periodic Solution#
This is a special case that might be needed, e.g. when simulating a section of a needle bearing with a 2D roughness profile, where we assume periodicity along the axis but non-periodic boundary conditions in the sliding direction.
Here, we numerically approximate the Green’s function by using the non-periodic solution and adding up periodic images in the direction where we want periodicity. For example, for periodicity in \(x\)-direction and non-periodicity in \(y\)-direction:
with the length of the domain in \(x\)-direction \(L_x\) and the number of summed periodic images in both positive and negative direction \(n\). For our purposes, we assume \(n = 10\) providing sufficient accuracy.
Fourier Transform Trick#
Since the convolution operation scales with \(N^2\) in real space (\(N\) being the number of grid points), we exploit the fact that in Fourier space, convolution is just a “simple” matrix multiplication. As transforming to Fourier space and back-transforming to real space can be done efficiently using Discrete Fourier Transform (DFT), this workflow scales only with \(N \cdot \text{log}(N)\) Stanley, Kato, J. Tribol. 119, 481 (1997).
Note that for the non-periodic and semi-periodic case, the technique of decoupling periodic images is required to use this Fourier Transform Trick. For that, we assume that the pressure outside the domain is zero. Refer e.g. to Hockney, Methods Comput. Phys. 9, 135 (1970) and Pastewka, Robbins, Appl. Phys. Lett. 108, 221601 (2016).
Usage#
YAML#
If you want to include elastic deformation in your model, you need to introduce the elastic field under properties. Therein, you can specify the material parameters
Young’s modulus \(E\) in Pa, defaults to 210 GPa
Poisson’s ratio \(v\) in - , defaults to 0.30
as well as solver-related underrelaxation parameter that dampens the elastic deformation response. If you experience issues after enabling elastic deformation, try making this value smaller.
alpha_underrelax, defaults to 1e-03
In the special case of semi-periodic boundary conditions, you can set the number of periodic images (evaluated in both directions, i.e. n_images=10 translates to 21 images being summed up to obtain the Green’s function).
n_images, defaults to 10
properties:
elastic:
E: 210e09
v: 0.3
alpha_underrelax: 1e-03
n_images: 10
Implementation Notes#
Periodicity is automatically derived from the density boundary conditions. If density is set periodic, elastic deformation will be calculated using the periodic Green’s function.
For non-periodic deformation, we define a reference pressure and reference deformation at the left side of the domain. This value is in both cases subtracted from all field values. Thereby, deformation is fixed to zero at the left boundary (This is a question of whether your given geometry should only define the initial state or if you want specific points to be fixed in the deformed state). See the
Parabolic Slider, Non-Periodicexample for clarification. You may want to change this fixation within theTopography.update()function withinGaPFlow/topography.py.Note that if you specify a non-periodic 1D problem in \(x\)-direction with \(N_y\)=1 and periodicity in \(y\)-direction, the solver will assume line contact for the calculation of the elastic deformation. That means that the non-periodic Green’s function is applied and \(L_y\) is set to the unit length of 1 m (This is relevant because the deformation is determined by the applied forces, and calculating the force requires specifying an area). As a consequence, this prevents accidental dependence of the result on \(L_y\).
Examples#
Journal Slider, Periodic#
This example shall highlight the deformation with periodic boundary conditions. Therefore we analyze a 1D journal slider and assume periodicity. Similar to the previous example, we assume a rather soft material with \(E = 5\) GPa to obtain significant deformation without requiring pressures in the GPa range.
sim = """
options:
output: data/journal
write_freq: 1000
silent: False
grid:
dx: 1.e-5
dy: 1.
Nx: 100
Ny: 1
xE: ['P', 'P', 'P']
xW: ['P', 'P', 'P']
yS: ['P', 'P', 'P']
yN: ['P', 'P', 'P']
geometry:
type: journal
CR: 1.e-2
eps: 0.7
U: 0.1
V: 0.
numerics:
CFL: 0.25
adaptive: 1
tol: 1e-8
dt: 1e-10
max_it: 40_000
properties:
shear: 0.0794
bulk: 0.
EOS: DH
P0: 101325.
rho0: 877.7007
T0: 323.15
C1: 3.5e10
C2: 1.23
elastic:
E: 5e09
v: 0.3
alpha_underrelax: 1e-04
"""
Now we initialize the problem
from GaPFlow.problem import Problem
myProblem = Problem.from_string(sim)
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- output : data/journal
- write_freq : 1000
- use_tstamp : True
- silent : False
- grid:
- Nx : 100
- dx : 1e-05
- Lx : 0.001
- Ny : 1
- dy : 1.0
- Ly : 1.0
- dim : 1
- bc_xE_P : [True, True, True]
- bc_xE_D : [False, False, False]
- bc_xE_N : [False, False, False]
- bc_xW_P : [True, True, True]
- bc_xW_D : [False, False, False]
- bc_xW_N : [False, False, False]
- bc_yS_P : [True, True, True]
- bc_yS_D : [False, False, False]
- bc_yS_N : [False, False, False]
- bc_yN_P : [True, True, True]
- bc_yN_D : [False, False, False]
- bc_yN_N : [False, False, False]
- geometry:
- U : 0.1
- V : 0.0
- type : journal
- flip : False
- CR : 0.01
- eps : 0.7
- numerics:
- tol : 1e-08
- max_it : 40000
- dt : 1e-10
- adaptive : True
- CFL : 0.25
- MC_order : 1
- properties:
- shear : 0.0794
- bulk : 0.0
- EOS : DH
- rho0 : 877.7007
- P0 : 101325.0
- C1 : 35000000000.0
- C2 : 1.23
- elastic:
- enabled : True
- E : 5000000000.0
- v : 0.3
- alpha_underrelax : 0.0001
- n_images : 10
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
Writing output into: data/2026-02-26_113452_journal
Let’s look at the initial geometry:
%matplotlib inline
myProblem.plot_topo()
Run the simulation and observe the result:
myProblem.run()
myProblem.plot_topo(show_defo=True, show_pressure=True)
-------------------------------------------------------------
Step Timestep Time CFL Residual
-------------------------------------------------------------
0 1.8984e-10 0.0000e+00 2.5000e-01 1.0000e+00
1000 1.8982e-10 1.8982e-07 2.5000e-01 3.4048e-04
2000 1.8981e-10 3.7964e-07 2.5000e-01 2.9804e-05
3000 1.8981e-10 5.6945e-07 2.5000e-01 7.9710e-06
4000 1.8981e-10 7.5927e-07 2.5000e-01 3.5214e-06
5000 1.8981e-10 9.4908e-07 2.5000e-01 1.5100e-05
6000 1.8981e-10 1.1389e-06 2.5000e-01 2.6030e-05
7000 1.8981e-10 1.3287e-06 2.5000e-01 3.5601e-05
8000 1.8981e-10 1.5185e-06 2.5000e-01 4.3432e-05
9000 1.8981e-10 1.7083e-06 2.5000e-01 4.9382e-05
10000 1.8981e-10 1.8981e-06 2.5000e-01 5.3484e-05
11000 1.8981e-10 2.0880e-06 2.5000e-01 5.5882e-05
12000 1.8981e-10 2.2778e-06 2.5000e-01 5.6785e-05
13000 1.8981e-10 2.4676e-06 2.5000e-01 5.6431e-05
14000 1.8981e-10 2.6574e-06 2.5000e-01 5.5059e-05
15000 1.8981e-10 2.8472e-06 2.5000e-01 5.2893e-05
16000 1.8981e-10 3.0370e-06 2.5000e-01 5.0137e-05
17000 1.8982e-10 3.2268e-06 2.5000e-01 4.6966e-05
18000 1.8982e-10 3.4167e-06 2.5000e-01 4.3530e-05
19000 1.8982e-10 3.6065e-06 2.5000e-01 3.9951e-05
20000 1.8982e-10 3.7963e-06 2.5000e-01 3.6331e-05
21000 1.8982e-10 3.9861e-06 2.5000e-01 3.2747e-05
22000 1.8982e-10 4.1759e-06 2.5000e-01 2.9261e-05
23000 1.8982e-10 4.3657e-06 2.5000e-01 2.5920e-05
24000 1.8982e-10 4.5556e-06 2.5000e-01 2.2757e-05
25000 1.8982e-10 4.7454e-06 2.5000e-01 1.9794e-05
26000 1.8982e-10 4.9352e-06 2.5000e-01 1.7047e-05
27000 1.8982e-10 5.1250e-06 2.5000e-01 1.4523e-05
28000 1.8982e-10 5.3148e-06 2.5000e-01 1.2222e-05
29000 1.8982e-10 5.5046e-06 2.5000e-01 1.0144e-05
30000 1.8982e-10 5.6945e-06 2.5000e-01 8.2807e-06
31000 1.8982e-10 5.8843e-06 2.5000e-01 6.6248e-06
32000 1.8982e-10 6.0741e-06 2.5000e-01 5.1654e-06
33000 1.8982e-10 6.2639e-06 2.5000e-01 3.8906e-06
34000 1.8982e-10 6.4537e-06 2.5000e-01 2.7878e-06
35000 1.8982e-10 6.6436e-06 2.5000e-01 1.8435e-06
36000 1.8982e-10 6.8334e-06 2.5000e-01 1.0446e-06
37000 1.8982e-10 7.0232e-06 2.5000e-01 3.7756e-07
37652 1.8982e-10 7.1470e-06 2.5000e-01 7.4811e-09
=================================
Total walltime : 0:02:41
(232.95 steps/s)
=================================
We can also generate an animation of the solution process (make sure you have the ipympl package installed):
myProblem.animate()
Note that, by using the Dowson-Higginson equation of state, we allow negative pressures in this example. For more physical behaviour, we should use an equation of state that considers cavitation and therefore does not output negative pressures. For that, see the non-periodic example below that uses Bayada-Chupin equation of state.
Parabolic Slider, Non-Periodic#
In this example, we want to illustrate how a parabolic slider deforms.
For this notebook, we will define this problem using a string variable as input. The structure is similar to the .yaml file.
This is a 1D problem with non-periodic boundary conditions, therefore the elastic deformation will be computed using the non-periodic Green’s function.
Note that we use a rather soft material with \(E = 5\) GPa to obtain significant deformation without requiring pressures in the GPa range.
sim = """
options:
output: data/pslider1d_elastic
write_freq: 1000
use_tstamp: True
grid:
Lx: 0.0762
Ly: 1.
Nx: 100
Ny: 1
xE: ['D', 'N', 'N']
xW: ['D', 'N', 'N']
yS: ['P', 'P', 'P']
yN: ['P', 'P', 'P']
xE_D: 850. # 101325
xW_D: 850. # 101325
geometry:
type: parabolic
hmin: 0.5e-5
hmax: 5.0e-4
U: 10.0
V: 0.
numerics:
adaptive: 1
CFL: 0.45
tol: 1e-8
dt: 1.e-10
max_it: 100_000
properties:
EOS: Bayada
rho0: 850.
shear: 0.039
bulk: 0.
cl: 1600.
cv: 352.
elastic:
E: 5e09
v: 0.3
alpha_underrelax: 1e-03
piezo:
name: Dukler
shearv: 3.9e-5
rhol: 850.
rhov: 0.019
"""
Now we initialize the problem
from GaPFlow.problem import Problem
myProblem = Problem.from_string(sim)
*************************************************************
* PROBLEM SETUP *
*************************************************************
- options:
- output : data/pslider1d_elastic
- write_freq : 1000
- use_tstamp : True
- silent : False
- grid:
- Nx : 100
- Lx : 0.0762
- dx : 0.0007620000000000001
- Ny : 1
- Ly : 1.0
- dy : 1.0
- dim : 1
- bc_xE_P : [False, False, False]
- bc_xE_D : [True, False, False]
- bc_xE_N : [False, True, True]
- bc_xW_P : [False, False, False]
- bc_xW_D : [True, False, False]
- bc_xW_N : [False, True, True]
- bc_xE_D_val : 850.0
- bc_xW_D_val : 850.0
- bc_yS_P : [True, True, True]
- bc_yS_D : [False, False, False]
- bc_yS_N : [False, False, False]
- bc_yN_P : [True, True, True]
- bc_yN_D : [False, False, False]
- bc_yN_N : [False, False, False]
- geometry:
- U : 10.0
- V : 0.0
- type : parabolic
- flip : False
- hmin : 5e-06
- hmax : 0.0005
- numerics:
- tol : 1e-08
- max_it : 100000
- dt : 1e-10
- adaptive : True
- CFL : 0.45
- MC_order : 1
- properties:
- shear : 0.039
- bulk : 0.0
- EOS : Bayada
- rho_l : 850.0
- rho_v : 0.019
- c_l : 1600.0
- c_v : 352.0
- rho0 : 850.0
- piezo:
- name : Dukler
- eta_v : 3.9e-05
- rho_l : 850.0
- rho_v : 0.019
- elastic:
- enabled : True
- E : 5000000000.0
- v : 0.3
- alpha_underrelax : 0.001
- n_images : 10
- gp:
- db:
- md:
*************************************************************
* PROBLEM SETUP COMPLETED *
*************************************************************
Writing output into: data/2026-02-26_113754_pslider1d_elastic
/usr/local/lib/python3.12/dist-packages/GaPFlow/topography.py:369: UserWarning: You specified a semi-periodic 1D problem.
For the calculation of elastic deformation, we assume a line contact with non-periodic boundary conditions in both directions.
For the calculation of the effective force F=p*A per cell, we assume a unit length of Ly = 1.
warnings.warn(
Let’s look at the initial geometry:
%matplotlib inline
myProblem.plot_topo()
Think about how there will be pressure building up on the left side of the slider and how the resulting deformation will change the shape.
Let’s run the simulation and plot the result:
myProblem.run()
myProblem.plot_topo(show_defo=True, show_pressure=True)
-------------------------------------------------------------
Step Timestep Time CFL Residual
-------------------------------------------------------------
0 1.9642e-07 0.0000e+00 4.5000e-01 1.0000e+00
1000 1.9194e-07 1.9058e-04 4.5000e-01 4.4222e-04
2000 1.9263e-07 3.8273e-04 4.5000e-01 4.4753e-05
3000 1.9259e-07 5.7527e-04 4.5000e-01 9.7512e-05
4000 1.9261e-07 7.6787e-04 4.5000e-01 7.1580e-05
5000 1.9261e-07 9.6048e-04 4.5000e-01 5.9350e-05
6000 1.9262e-07 1.1531e-03 4.5000e-01 5.4108e-05
7000 1.9262e-07 1.3457e-03 4.5000e-01 5.0155e-05
8000 1.9261e-07 1.5383e-03 4.5000e-01 4.6667e-05
9000 1.9261e-07 1.7309e-03 4.5000e-01 4.3203e-05
10000 1.9261e-07 1.9236e-03 4.5000e-01 3.9839e-05
11000 1.9262e-07 2.1162e-03 4.5000e-01 3.6702e-05
12000 1.9262e-07 2.3088e-03 4.5000e-01 3.3563e-05
13000 1.9261e-07 2.5014e-03 4.5000e-01 3.0344e-05
14000 1.9261e-07 2.6940e-03 4.5000e-01 2.7146e-05
15000 1.9261e-07 2.8866e-03 4.5000e-01 2.3454e-05
16000 1.9261e-07 3.0793e-03 4.5000e-01 2.1021e-05
17000 1.9261e-07 3.2719e-03 4.5000e-01 1.7340e-05
18000 1.9261e-07 3.4645e-03 4.5000e-01 1.6156e-05
19000 1.9261e-07 3.6571e-03 4.5000e-01 1.2202e-05
20000 1.9261e-07 3.8497e-03 4.5000e-01 1.1255e-05
21000 1.9261e-07 4.0423e-03 4.5000e-01 8.1897e-06
22000 1.9261e-07 4.2349e-03 4.5000e-01 6.8245e-06
23000 1.9261e-07 4.4276e-03 4.5000e-01 4.5701e-06
24000 1.9261e-07 4.6202e-03 4.5000e-01 4.2918e-06
25000 1.9261e-07 4.8128e-03 4.5000e-01 2.4779e-06
26000 1.9261e-07 5.0054e-03 4.5000e-01 2.0989e-06
27000 1.9261e-07 5.1980e-03 4.5000e-01 7.4475e-07
28000 1.9261e-07 5.3906e-03 4.5000e-01 1.6563e-07
28109 1.9261e-07 5.4116e-03 4.5000e-01 2.4696e-09
=================================
Total walltime : 0:02:20
(200.64 steps/s)
=================================
Note here, as discussed in Implementation Notes, that at the left boundary, the deformed shape is fixed to the initial shape (meaning the effective deformation is set to zero here). With no fixing, the deformation would be positive for the whole domain and the whole slider would be pressed upwards.